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ABSTRACT 


The modern aerospace controls engineer is provided with a variety of powerful 
tools to aid in the design and testing of digital flight control systems. The current fis¬ 
cal environment requires extensive validation of all aerospace based systems through 
simulation and hardware-in-the-loop testing prior to implementation. This work ex¬ 
plores the design and evaluation of an Automatic Flight Control System (AFCS) for 
the Bluebird Unmanned Aerial Vehicle (UAV). Software tools such as Matlab and 
MATRlX;i' are used to evaluate the dynamic stability of the aircraft model and Linear 
Quadratic Gaussian algorithms are used to obtain the appropriate controller. Graph¬ 
ical design applications such as SiMULiNK and SystemBuild are then used to build 
a visual block diagram model of the aircraft dynamics and link it with the designed 
controller. Using this model, the control system response to commanded input.s and 
external disturbances was evaluated. 
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I. INTRODUCTION 


A. BACKGROUND FOR UAV RESEARCH 

The modern battlefield has increasingly progressed towards the use of automated 
systems and remotely controlled devices to perform a variety of missions. From 
surveillance to weapons delivery and bomb damage assessment, the human operator 
is being remo^’ed from the direct danger of a hoslile environment and placed in 
a position of evaluating data received via. RF or fiber optic link. The direct and 
obvious benefits of such an arrangement cire the reduced risk to the operator and the 
reduced cost of the unmanned sensor platform as compared to traditional manned 
platforms. The state-of-th^art technology in unmanned aerial vehicle development 
has demonstrated the capability of flight out to ranges of .500 nm and endurances 
exceeding 24 hours. Combined with the ability to carry a variety of sensor suites, 
these platforms represent the future in airborne data acquisition for both military 
and civilian applications. 

In support of these technological developments the Unmanned Air Vehicle Flight 
Research Lab (UAV FRL) at the Naval Postgraduate School has been investigating 
several unmanned aerial vehicles as technology demonstrators. The .AROD U.-W is a 
vertical take-off and landing platform. Vertical flight is accomplished with a powerful 
ducted fan producing enough thrust to lift the aircraft. Current proposals have 
the AHOD working in conjunction with unmanned surface and subsurface vehicles 
providing additional remote sensing capabilities and data link services between the 
operator and the surface vehicles. The AROD is an inherently unstable platform 
and is subject to gyroscopic coupling and torque effects during the production of 



TABLE 1.1; PHYSICAL CHARACTERISTICS OF Bluebird 


Weight 

55 lbs 

Average Wing Chord, c 

1.802 / 

Wing Span, b 

12.42 / 

Planform Surface Area, S 

22.380 p 

Engine Power 

4.0 HP 

Mass Moment of Inertia, 

10.0 slug - p 

Mass Moment of Inertia, ly 

16.12 slug - p 

Mass Moment of Inertia, 1^ 

7.97 slug - p 


lifting thrust. Extensive modeling and simulation of this vehicle was previously 
accomplished by Sivashankar and Moats [Ref. 1, 2] in separate work at the UAV FRL. 
This work validated the dynamically unstable nature of the AROD and provided the 
motivation for the second UAV project currently under development at the UAV 
FRL. 

The Bluebird aircraft was acquired as a test bed for guidance and navigation 
systems. It is similar in appearance to a scaled down Cessna 172. It’s physical 
characteristics are given in Table 1.1. Its conventional high-wing configuration makes 
for a stable aircraft. This provides the ideal platform for testing guidance, navigation, 
and control software and hardware before installation on the AROD. As with the 
AROD, the Bluebird has been extensively modeled [Ref. 3, 4], the results of which 
will be covered in Chapter III. 

B. OVERVIEW 

This thesis fulfills a twofold purpose. First, to provide for the autonomous con¬ 
trol of the Bluebird U.AV, a controller is designed based on Linear Quadratic Regula¬ 
tor Theory and using the ’Iqr’ and ’regulator’ functions of Matlab and MATRIX^- 
This design will allow the remote operator to control the vehicle’s altitude, true air- 
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speed, and heading, while limiting the response to commanded inputs to within the 
vehicles dynamic operating envelope. This stable control provides the capability to 
test a variety of avionics systems through a range of dynamic maneuvers that would 
not be possible in tethered flight. It also provides for a more stable control than in the 
case of direct RF control by a human operator. Second, this work will provide a link 
from the courses in classical and modern control theory at the Naval Postgraduate 
School to the implementation of these concepts using state-of- the -art software tools 
such as MATL.4B and MATRIX;^. The ultimate goal is the integration of the model¬ 
ing of the airborne platform and sensors, controller design, and hardware--in-the-loop 
testing of the design on the chosen platform. 

These objectives wcire achieved in a multi-step process described in this thesis. 
This description begins with a summary'of the development of the nonlinear equations 
of motion of a rigid body in space that is subjected to external forces and moments. 
The formulation of these equations has been the subject of much study at the UAV 
FRL [Ref. 3, 4]. For this reason only the significant results will be examined as an 
aid in understanding the development of the aircraft models (Chapter 11). 

Following this step, the equations of motion are encoded to form the core of a 
high fidelity nonlinear block diagram model of the aircraft dynamics in SlMULINK and 
SystemBuild. These computer codes have been previously developed and validated 
independently [Ref. 4, 5] in the .m file format of Matlab-Simulink and as C- 
code. lb model disturbances induced during flight through a moving air mass and 
to calculate the aircraft true airspeed, this work modified these computer codes to 
include wind inputs in the inertial x, y, and coordinate directions. 



To better determine control system requirements an open-loop analysis of the 
aircraft model was done as follows: 

• The nonlinear model was trimmed about a nominal operating point around 
which the dynamic response to small perturbations could be analyzed. 

• The model was then linearized around the trim point to obtain a linear model. 

• The eigenvalues of the linear plant were determined and the natural frequency 
and damping of the different modes were analyzed. 

. From this data controller requirements were established, determining desired band- 
withs for response to command and control inputs versus the actual open-loop plant 
responses (Chapter III). 

These requirements provide the basis evaluating the feedback controller obtained 
using the linear quadratic regulator algorithms of MatLAB and MATRIX a;. The 
controller design was based upon linear quadratic regulator theory. To allow for a 
better understanding of the algorithms used to calculate the controller, the main 
points of this theory are reviewed. The controller design proceeds using the following 

• The control synthesis model is developed. In this model the states to be con¬ 
trolled and actuators to accomplish this control are included. 

• The control gains are calculated using the appropriate MaTLAB and M.ATRIXy 
algorithms. In these calculations a cost function, which includes weighting 
factors, is used to modify the energy penalty incurred in responding to the 
various control and command inputs. 



• The time and frequency domain response of the closed-loop controller and 
plant are evaluated. If the step response or bandwith of the system does not 
meet the design requirements, the weighting factor is changed for the particular 
command or control input and the process* repeated. 

The entire design process may take several iterations and may also involve the ma¬ 
nipulations of zeros added to the synthesis model. In theory, as the weights on 
the controls go to zero, the closed loop poles will go to the open-loop zeros that 
were added to the model. This ensures that the closed-loop poles have the desired 
damping ratio and natural frequency to obtain the desired system response (Chapter 
I\'). 

Once the controller meets design requirements, it is implemented on the nonlin¬ 
ear plant. This implementation is accomplished in four steps: 

• Linear controller with the linear plant, initial conditions equal to zero. 

• Differential controller with the linear plant, initial conditions equal to zero. 

• Differential controller w’ith the nonlinear plant, initial conditions equal to the 
trim states. 

• Differential controller and nonlinear plant are discretized to run on a digital 
computer. 

At the completion of each stage of this integration process the closed-loop con¬ 
troller/plant model is tested to ensure the appropriate response to a variety of com¬ 
manded inputs and external disturbances (Chapter V). 



II. BACKGROUND: EOM DEVELOPMENT 


The development of the equations of motion for a generic aircraft with six degree 
freedom of motion has been the subject of much study at the U AV FRL . High fidelity 
nonlinear models of both the AROD and Bluebird have been developed and verified 
through extensive simulation [Ref. 3, 4], To rederive these results would be an 
unnecessary repetition of effort. However, to provide a sound basis for understanding 
the models developed in this thesis, a summary of the significant results of these 
previous works will be provided. 

A. FORCE EQUATION 

The derivation of equations of motion for a general six degree of freedom airplane 
model can be divided into two parts. The first part is simply the determination of the 
equations of motion for any rigid body in space. It is dependent only on the linear 
and angular momenta of the body. The equations of linear motion for the aircraft 
center of gravity are derived by a direct application of Newton’s Law, F — ma. The 
final result for the forces acting on the aircraft body in a body referenced coordinate 
system is, 

= m {—^vbo + X ^i>Bo) 

= rn -^^vbo + in ^ujb x '^vbo- (2.1) 

• represents the forces on the aircraft in the body coordinate system. 
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• represents the derivative with respect to time in the body coordinate 
sytem of the velocities of the aircraft body origin. 

• are the velocities of the body origin resolved in the body coordinate 
system. 

• ^lob are the body angular rates resolved in the body coordinate system. 

B. MOMENT EQUATION 

The equations for angular acceleration are derived from Euler’s Law for the 
conservation of angular momentum, 

i = iV, (2.2) 

The final result, shown in Equation 2.3, expresses the total external moment 
applied to the aircraft center of gravity and is given as: 

^Nbo = Ib^^b + X Ub^i^b + (2.3) 


• represents the moments acting on the aircraft body origin resolved in 
the bodj' coordinate system. 

• Ib Is the inertia tensor for the aircraft. 

• Ir is the inertia tensor for any significant rotating object located at the center 
of gravity of the aircraft, such as a propeller or a turbine. It should be noted 
that the term Ir^ijJr can be disregarded if it is insignificant compared to Ib 
and ^uJb [Ref. 6]. 




C. EXTERNAL FORCES AND MOMENTS 


The second part is the calculation of aerodynamic, gravitational, and thrust 
forces and moments on the airplane. The total forces and moments acting on the 
aircraft are determined by summing these components and can be expressed as: 

1 — [ + ^FpROP + AERO 

J - [ ^Nprop + ^Naero 

The aerodynamic forces are specific to the aircraft platform and are calculated 
using nondimcnsional stability and control derivatives. The stability and control 
derivatives are obtained by approximating the aerodynamic forces and moments act¬ 
ing on the aircraft using a first order Taylor Series expansion about a given flight 
trim point [Ref. 7]. Since these derivatives are generally computed in the wind ref¬ 
erenced coordinate system it is necessary to resolve the resulting expression for the 
aerodynamic forces and moments in a body referenced coordinate system. The final 
result for the aerodynamic forces is given by [Ref. 8]; 


(2.4) 


{Cfo + -^M' X -f -^M'i + ^A} 


• q is the dynamic pressure (0.5pF*). 

• 5 = diag{5, S, S, Sb, Sc, Sb}. 

• A = [SeUvatoT^ brudder, baiieron] > is the vector of Control inputs. 

• is the rotation matrix from wind to body coordinates. 

• M' is the scaling matrix given by diag{l/Rj, 1/Vr, 1 /Vr, b/2VT, c/2Vt, 6/2Vr}' 

• Cfo is the vector of steady state coefficients representing trimmed flight. 
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Ip is the matrix of stability and control derivatives and is given by, 


Cl^- 

Cl, 

Cl. 

CLp 

Cl, 

Cl. 

Cyu 

Cy, 

Cy. 

Cyp 

Cy, 

Cy. 

Cd, 

Cd, 

Cd. 

Cd. 

Cd, 

Cd. 

Ci, 

Cl, 

Cl. 

Cip 

Cl, 

Cl. 

Cmu 

Cm, 

Cm. 

Cmp 

Cm, 

Cm. 

Cnu 

Cn, 

Cn. 

Cn. 

Cn, 

Cn. 


( 2 . 6 ) 


Ip is very similar to |p, except that only the a and ^ terms are the only 
nonzero terms. 


• X is the state vector composed of body inertial velocities and rates. 

The calculation of the forces and moments due to gravitational and propulsive 
efFe( Is are carried out in a more straight forward manner. The final results are given 
bt-. 

[ ° 1 

^ Foray = 0 (2.7) 

L '<^9 J 

^FpRop — Ty (2.8) 

[Tz\ 

and 

r T> ] 

^NpRop = Tm , (2.9) 

L Tn J 

where the Ti terms represent the forces and moments due to the generation of thrust. 

D. STATE SPACE REPRESENTATION 


From the expressions derived for the forces and moments acting on the aircraft 
it is possible to derive a state space representation of the equations of motion. We 



can write equations 2.1 and 2.3 in the following manner, 

[ 1 _ r m jl^VBO + m x ‘^vbo) 1 

[ ^N \ - [ Ib’^^b + X + Ir^Ur) J 

By rearranging this expression and normalizing by 1/m and /g' \ 

suiting equation in state space form, 

£\ 

dt 

Substituting the expressions for the aerodynamic, gravitational, and thrust 
forces and moments into equation 2.4 and then substituting this into equation 2.11 
the complete state space representation of the nonlinear equations of motion becomes, 

' '' \ -®u)bX 


f ^VBO ] 

[ X ^VBO 

+ ^ 1 

[ J = 

[ -/b' Scub X {Jb^ujb + Ir^o^h 

) + 


( 2 . 10 ) 
? obtain the re- 


(2.11) 


dt 


= x- 


0 


Mr^fyTqS^M' J ^ 

r ^FpRop 

I ^NpRop 


X (^Ib^^b + Ir^<^r)) j 

Foray 1 


vbo 


h MT^ ■ 


j 8t + iTqSiCp, + ^A) j 


X' h ■ 

Ml = I 


Mj-^^,TqSS§ 

" 1 

) ®/b j 

%R 0 1 


0 ^,R J 

To complete the equations of motion, the following t 
e used to calculate the body positions and Euler angles, 

^Pbo = %R^vbo, 


3 differential equations 
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• gi? is the rotation matrix from the body to the inertial coordinate system. 

• i'tA) is the rotation matrix that takes the body angular rates to the Euler 
angular rates. 
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III. OPEN-LOOP ANALYSIS 


In order to design a controller which stabilizes the feedback system for the 
Bluebird UAV it is necessary to analyze the open-loop characteristics of the aircraft 
model. Before this analysis takes place, the differential equations of motion must be 
modeled using such tools as SiMULINK or SystemBuild or encoded in a form that 
graphical software applications such as these can use. Work by Halberg [Ref. 4] 
and Byerly [Ref. 5] have developed such codes in the .m file format of MaTLAB and 
as a C — code file in the User Code Block format of MATRIX.V. To account for 
the motion of the aerodynamic body through a moving airrnass, the C — code was 
modified to allow for the input of wind disturbances in the three inertial coordinate 
directions, The .tn file allows for airniass disturbances and no modifications were 
necesSMy, TbC: codes are shown in .Appendix A. Using these codes as the core, both 
SlMfLlNK and SystemBuild block diagram models were developed to represent the 
dynamic aircraft model. Using data for the Cessna 172 aircraft from Roskam [Ref. 
6] both the models were validated by comparing eigenvalues for the open-loop plant. 
Using these models the open-loop analysis proceeded. This analysis provided data 
on the damping and frequencies of the different aircraft modes. Additionally, the 
time history of the aircraft states was obtained and plotted to provide a visual clue 
to the open- loop aircraft performance. Using this data, control requirements were 
derived and controller design accomplished. 
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A. PLANT MODELING AND VALIDATION 


The computer codes modeling the aircraft equations of motion have been pre¬ 
viously developed. The validation of these codes using SiMULINK or SystemBuild 
requires the sreneration of anr>ronria.te block diafframs. Fieiire .1.1 and sho^vs the 
structu e block 



Figure 3.1: BlueBird Open-loop Plant Model - SiMULINK [Ref. 4] 

named EOM/wind, Figure 3.2, shows the actual implementation of the nonlinear 
differential equations of motion. The M.^TLAB Function block named ’state_dcriv’ 
is the calling block for the .m file that calculates the derivative states of the body 
velocities and angular rates. 

The MATRIX;^ implementation of the equations of motion was less complex 
since the calculation of the position and Euler angle derivative states was internal to 
the C - code routine that is called in the SuperBlock USR003. Figure 3.3 details the 
SystemBuild open-loop model. Note that the only external routine block required is 
for the calculation of the vehicle true airspeed. 

To determine the vehicle open-loop characteristics and performance, so that a 
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Figure 3.2: SiMULIN'K Implementation of the EOM, [Ref. 4] 

linear controller can be designed, both models were trimmed using the appropriate 
Matlab and MATRIX;f functions and then linearized about the trim state. The 
open-loop plant is shown in Appendix B. According to small perturbation theory, it 
can then be assumed that the resulting aircraft plant will respond linearly to small 
magnitude disturbances about the trim state. 

To validate the model dynamics, stability and control derivatives from Roskam 
[Ref. 6| were substituted into the coded equations of motion. The models were then 
trimmed and linearized. The eigenvalues of the resulting linear state space model were 
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Figure 3.3: SystemBuild Open-Loop Plant Model 


then determined and compared with the values from Roskam. The results are shown 
in Table 3.1. The data in Table 3.1 indicates that the models are highly accurate. 
The only significant difference occurs in the real part of the Phugiod mode poles. 
However, when the magnitude of these eigenvalues are compared, the difference is 
less than 0.2%. 

With the models validated, the open-loop analysis of the Bluebird plant model 
followed the identical procedure. After substituting the aerodynamic stability and 
control derivatives in the encoded equations of motion, the model was then trimmed 



1 Mode 

Roskam 

Matlab 

M.ATRIXx ' 

Longitudinal 




Short Period 

-4.130 + 4.390i 

-4.1303 -1- 4.3895; 

-4.1303 -f 4.3895; 


-4.130 - 4,390i 

-4.1303 - 4.3895; 

-4.1303 - 4.3895; 

Phugoid 

-.02092 + .1797j 

-.0135 + .1801; 

-.0135-h .1801; 


-.02092 - .1797i 

-.0135 -.1801; 

-.0135 - .1801; 

Lateral-Directional 




Dutch Roll 

-.6858 + 3.306i 

-.6930 -h 3.3077; 

-.69.30 + 3.3077; 


-.6858 - 3.306; 

-.6930 - 3.3077; 

-.6930 - 3.3077; 

Roll 

-12.43 

-12.4343 

-12.4343 

Spiral 

-.01095 

-.0109 

-.0109 


TABLE 3.1; CESSNA 172 EIGENVALUES 
and linearized. The results were identical for the SiMULiNK and SystemBuild models. 
These results are summarized in Table 3.2 along with the associated damping ratio 
and frequency of each mode. 


! Mode 

Eigenvalue 

Damping, ( Frequency(rad/sec), jj \ 

Longitudinal 




Short Period 

-4.3290 4- 3.9939; 

.7.35 

5.8899 


-4.3290 - 3.9939; 

.735 

5.8899 

Phugoid 

-.0171 + .4970; 

.0344 

.4973 


-.0171 - .4970; 

.0344 

.4973 

Lateral-Directional 




Dutch Roll 

-.2665 + 2.3861; 

.1110 

2.4009 


-.2665 -2.3861; 

.1110 

2.4009 

Roll 

-4.5722 

1.0 

4.5722 

Spiral 

.0384 

-1.0 

.0384 


TABLE 3.2: Bluebird EIGENVALUES 


To provide further data for the open-loop analysis, the model simulations were 
run over a period of 300 secs and time history plots of the output states were obtained. 


(Appendix B) 




B. CONTROLLER DESIGN REQUIREMENTS 

The information provided by the open-loop analysis lead to several conclusions 
about the uncontrolled aircraft dynamics: 

• The short period is critically damped with stable left plane poles. 

• The phugoid mode is extremely underdamped. To improve the altitude and 
airspeed tracking performance of the vehicle, the damping of this mode must 
be increased. 

• The dutch roll mode damping is inadequate for stable heading control. 

• The unstable spiral mode has a destabilizing effect on performance. Even 
though this mode is extremely slow, to adequately control aircraft heading 
this eigenvalue must be moved to the left half plane. 

Additional considerations in establishing controller requirements for the Blue¬ 
bird I.IAV were: 

• The Bluebird test vehicle is designed for use as an avionics test bed. The plat¬ 
form flight regime should thus include only relativeh’ benign manuevers. This 
requirement leads to the implementation of limiters on both body coordinate 
X— and z— accelerations and bank angle, 4>- 

• The Futaba SB-34 servomotors that serve as control actuators have been shown 
to respond accurately up to 4 rad/sec [Ref. 2]. 

The preceding factors lead to the formulation of the following control system 
requirements: 

• Feedback system must be stable. 
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Minimum phugoid mode damping of 0.50. 


Yaw rate commanded continuously to zero. 

Control bandwidths less than or equal to 4radlsec to conform with the servo¬ 
motor limitations. 

Maximum overshoot of 10% to step commands in airspeed, altitude, and head¬ 
ing. 

Rise lime to step commands in airspeed, altitude, and heading of approximately 
4sec. Rise time is defined as the time to go from 10% to 90% of the commanded 
input. 

Maximum angle of bank of 30 degrees. 

Accelerations in the x— and z— body coordinates limited to less than I.O5. 



IV. LQR THEORY AND CONTROLLER 
DESIGN 


A. OPTIMAL CONTROL 

In designing control systems we are often concerned with choosing the control 
u(i) such that a given performance index is minimized. Such an index is given by: 

J = y {x'Qx + u'Ru) dx (4.1) 

Now suppose a system is described by the following equations; 

X = A X B u 

y = 0. (4.2) 

assuming that (.4. B) is stabilizable, (Q.A) is detectable, Q>0, and i? > 0. 

From theory [Ref. 9], given the performance index shown in Equation 4.1 and 
the system gi\’en by Equations 4.2, there exists a minimizing control, 

u = K X (4.3) 

which stabilizes {A — BK), w'here: 

K=R-^B'P (4.4) 

In Equation 4.4, P is a unique, positive definite solution to the Algebraic Ricatti 
Equation shown in Equation 4.5. 


A'P + PA- PBR-^B'P + Q = 0 




By solving this equation for the matrix P, for the chosen weighting matrices Q and 
R, the optimal gain K can be determined using equation 4.4. 

The entire control design process reduces to an iterative procedure that follows 
these basic steps: 

• Choose weighting matrices Q and R. 

• Calculate the solution to the Ricatti Equation, P, equation 4.5. 

• Calculate the optimal gain, K, equation 4.4. 

• Examine the closed-loop system for the desired frequency and step responses. 

• If the bandwith or time domain characteristics are not as desired, choose new 
Q ox R and begin the design procedure again. 

In this manner the designer is able to iterate through the family of stabilizing con¬ 
trollers until obtaining the desired frequency and time domain characteristics. 

B. CONTROLLER DESIGN AND ANALYSIS 

The control system was designed assuming zero coupling between the lateral and 
longitudinal states. Examination of the open-loop plant in .Appendix B shows the 
absence of cross-coupling terms in the state matrix, A, and minima! cross coupling 
in the control matrix, B. This analysis supports the choice of developing separate 
controllers for the lateral and longitudinal states. 

The previous section showed that the solution to the Ricatti Equation leads to 
the calculation of the optimal control based upon the choice of weighting matrices. To 
aid the modern controls designers, software packages such as MatlaB and MATRIX,v 
have provided algorithms that automate this process. These routines take as inputs 
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the plant and weighting matrices and output the optimal gain matrix, solution to 
the Ricatti equation, and the eigenvalues of the closed-loop plant. 

In the first step of utilizing these algorithms the open-loop plant is separated 
into longitudinal and lateral parts. The longitudinal plant is composed of the xi, xu, 
q, 6, and z states. The lateral plant is composed of the v, p, r, 4>, and 4' .states. The 
control synthesis models for the longitudinal and lateral plants were then developed. 
In these models the states to be controlled and the actuators to achieve this control 
were chosen. The longitudinal synthesis model is shown in Figure 4.1. The states 
to be controlled in the steady state were the inertial velocity, u, and the altitude, z. 
The actuators used to achieve this control were the elevator and throttle. The lateral 



Figure 4.1: SiMULiNK Longitudinal Control Synthesis Model 

state to be controlled in the steady state was heading, xj). The actuators recruited to 
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control V' were the ailerons. Additionally, to increase lateral stability, a yaw damper 
was implemented by feeding back yaw rate through a gain block directly to rudder 
input. The choice of gain is arbitrary and the final value was obtained via a trial 
and error process by examining the step response of the lateral states over a range of 
gains. The lateral synthesis model is shown in Figure 4.2. From an examination of 
the closed-loop from rudder input to yaw rate it can be shown that the time response 
of the yaw rate is a decreasing exponential function of the form; 

^ = Cl (4.6) 

where the constants, ci and C2, are elements of the open loop A and B matrices. 
Therefore this implementation of the yaw damper continuously commands yaw rate 
to zero. 



Note the addition of zeros to all the states being controlled. In theory, as the 
fights on the controls go to zero, the closed-loop poles will go to the open-loop 





zeros of the synthesis model. This ensures that the closed-loop poles have the desired 
damping ratio and natural frequency to obtain the desired system response. 

Once the control synthesis models were built, they were linearized using the ’lin- 
mod’ function of Matlab. -Along with the choices for the weighting matrices, these 
linear plant models provided the inputs to the ‘lqr‘ function of Matlab. This routine 
calculates and outputs the optimal gain along with the solution to the Ricatti equa¬ 
tion and the closed-loop eigenvalues. An .m file was written to compute the control 
gains for both the longitudinal and lateral controllers and is detailed in Appendix C. 
Utilizing a design technique presented in [Ref. 8], the optimal gains were multiplied 
by a factor of two. This had the effect of surpressing the magnitude of the frequency 
response at the resonant frequency, while slightly increasing the bandwidth of the 
response. The exact mechanics of this effect are unclear and are recommended for 
further study. This program also computes and plots the frequency response for the 
command and control loops of the controller as well as the time domain response to 
step command inputs. These plots are also detailed in Appendix C. 

Table 4.1 lists the closed-loop eigenvalues. .All modes are stable and mini¬ 
mum longitudinal damping is only slightly less than the desired value of 0.5. Fig¬ 
ures C.3, C.o, C.8 in Appendix C detail the respori.S)’#. to step altitude, velocity, and 
heading commands. Both the altitude and velocity responses, Figures C.3, C.5, easily 
meet the 4.0 seconds rise time and show no (Avrshoot to the step commands. The 
heading response overshoot, Figure C.8, is well under the established requirement at 
less than 1.0 

Finally, Figures C.l, C.2, C.7 of Appendix C indicate that the controller is 
keeping the actuator responses at or below the desired 4 rad/sec bandwidth, with 
minimal magnitude disturbances at the resonant frequencies. 
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1 Mode ^ Eigenvalue | 

Damping, ^ 

Frequency(rad/sec), lo | 

Longitudinal 





-4.3144 +4.0187i 

.7317 

5.8961 1 


-4.3144 - 4.0187j 

.7317 

5.8961 1 


-1.0933 

1.0 

1.0933 I 


-.4846 -h .8640y 

.4892 

.9906 1 


-.4846 - .8640; 

.4892 

.9906 


-.2692 -t- .2580; 

.7219 

.3728 


-.2692 - .2580; 

.7219 

.3728 

Lateral-Directional 





-1.2874 + 1,8991; 

.5611 

2.29:3 ! 


-1.2874 - 1.8991; 

.5611 

2.2943 


-.4347 + .5657; 

.6093 

.7134 


-.4347 - .5657; 

.6093 

.7134 


-4.5604 

1.0 

4.604 1 


-.65.56 

1.0 

.6556 1 


TABLE 4.1; Bluebird CLOSED-LOOP EIGEN\^LUES 
MATRIXjr also has the appropriate algorithms to calculate an optimal controller 
using the solution to the Ricatti Equation. The above results based on Matlab and 
SlMfLlNK were verified by Byerly [Ref. 5] using the the MATRIX^- codes with 
SystemBuild block diagrams. These results will not be repeated here. However, to 
allow for hardware-in-the-loop testing with the Futaba actuators, the controller will 
be implemented in SystemBuild as well as SiMULINK. This work is discussed next. 
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V. CONTROLLER IMPLEMENTATION AND 
TESTING 


Tbe previous chapter detailed the design and analysis of the linear controller. 
The next step is to implement the controller with the aircraft block diagram model 
and subject the feedback system to external commands and disturbances. The soft¬ 
ware packages of Matlab and MATRIX^ each offer graphical design tools, StviULLN'K 
and SystemBuild respectively, in which to implement the controller. The Matlab 
/ S 1.MULINK combination offers the advantage of familiarity. It is the standard engi¬ 
neering software for the .4eronautics and Astronautics Department at the Naval Post¬ 
graduate School. The nonlinear equations of motion for a vehicle moving through 
three-dimensional space have been well developed in the .m file format of Matlab 
[Ref. 3, 4] and SlMULlNK allows for the easy implementation of these user defined fun- 
tions into the block diagrams representing the aircraft plant. The major disadvantage 
is the inability to generate autocode. The MATRIX a- / SystemBuild combination 
allows the user to automatically generate C-code from the system block diagram. 
This computer code can then be run on any number of digital computers that are 
currently available. The disadvantage of using this software is that it is not as widely 
used and therefore not as familiar as Matlab. 

As noted above, each software package and its associated graphical design ap¬ 
plication have their advantages and disadvantages. For this reason the controller was 
first implemented using SlML’LINK. This implementation was accomplished in four 

• Linear controller with the linear plant, initial conditions equal to zero. 
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Differential controller with the linear plant, initial conditions equal to zero. 


# Differentia] controller with the nonlinear plant, intial conditions equal to the 
trim states. 

• Discretize the differenticJ controller and nonlinear plant to be compatible with 
a digital control computer. 

Once the aircraft/controller model was verified, the complementary model was de¬ 
veloped using SystemBuild and verified once again. Then the SystemBuild model 
was discretized and autocode generated. 

A. CONTROLLER STRUCTURE 

The optimal control gain, K, was obtained using the ’Iqr’ algorithm of Matlab. 
This algorithm used a linearized model of the synthesis model as part of its input 
variables. These synthesis models used integrators to develop the integral error be¬ 
tween the actual state and the commanded state. Thus the control gain can be 
separated into proportional and integral gain matrices and the controller must have 
both proportional and integral error portions. The proportional part assumes all the 
state variables are measurable. These variables are then used for state feedback . 
The integral portion consists of integrators on the error signal generated between 
the actual state and and the commanded state. The state space equations of the 
controller take the following form; 

Xc = Bc{y-yc) (5.1) 

u = Kix^ + K^x 

where, 
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• y = J vector of actual outputs. 

• = I j , is the vector of commanded outputs. 

• K = [/VpA'.] 

It is from this form that the SlMULINK implementation of the controller was devel¬ 
oped. The details of this development are covered in the following sections. 

B. LINEAR CONTROLLER, LINEAR PLANT MODEL 

The linear controller is shown in Figure .5.1. The controller takes the optimal 



Figure 5.1: SlMULINK Linear Control Model 


gain matrix calculated using the ’Iqr' of MaTLAB and breaks it into a matrix of 
proportional gains corresponding to the states and a matrix of integral gains corre¬ 
sponding to the integral error states. The control outputs were the commands for the 
longitudinal controls (elevator and throttle), the lateral control (ailerons), and the 
yaw damper (rudder). The complete closed-loop linear control/linear plant model 



is shown in Figure 5.2. The block named Long.Lat.Ctrl contains the linear control 
detailed in Figure 5.1. 



Figure 5.2: Simulink Closed-Loop Linear Control Model 

C. DELTA CONTROLLER, LINEAR PLANT MODEL 

The linear controller model of the previous section was designed for the linearized 
model of the aircraft dynamics. This linear plant model was obtained by trimming the 
nonlinear model about some trim flight state and then linearizing the model at that 
state. This entire process is performed under the assumption of small perturbations 
about the trim slate. Thus the linear controller on the linear plant is valid only 
for minor deviations about the trim position. It is then necessary to obtain small 
perturbations of the full states. This is achieved by differentiating the states prior 
to multiplication by the proportional gains and then integrating at the controller 







output. 

The complete delta controller is implemented as shown in Figure 5.3. The 
integrator has been removed from the integral error states and a differentiator has 
been added to the full states prior to the proportional gains. Once these signals have 
been summed the result is fed through an integrator, whose initial conditions equal 
the trim states. The resulting outputs are the control commands. 



Figure 5.3: Delta Controller Model 

D. DELTA CONTROLLER, NONLINEAR PLANT MODEL 

The linear controller used a combination of throttle and elevator to control 
inertial velocity, u. This control is extremely accurate for the case of an aerial 
vehicle moving through a still air mass, as the vehicle true airspeed (TAS) and inertial 
velocity (ground speed) are identical. However, for motion in a moving air mass it 
is desired to control the vehicle TAS. This is necessary to prevent stall. TAS is 
computed in two steps: 
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• Wind vectors in the three inertial coordinate axis are used as inputs to the 
nonlinear equations of motion. This includes the effects of air mass velocities 
in the calculations of the vehicle states. 

• The inertial velocities are then summed with the wind velocity components and 
the magnitude of their vector sum is determined by 

TAS = + + (.5.2) 

where, Wy, and are the wind components in the three inertial directions. 

This calculation is shown in the ’Fen’ block of Figure 3.2 and in the block TAS 
C.4LC’ of Figure 3.3. The inertial velocity, u. is then replaced by TAS to determine 
the velocity error. The actual controller implementation remains the same as shown 
in Figure 5.3. The closed-loop plant takes the form as shown in Figure 5.4, where the 
block ’EOM/wind' contains the implementation of the nonlinear equations of motion 
shown in Figure 3.2. 

E. DISCRETE CONTROLLER IN Simulink AND System- 
Build 

From the outset, the controller designed in this thesis was intended to operate 
on a digital flight control computer. For this reason the delta analog controller from 
the preceeding section must be discretized. This process involves implementing the 
controller in the discrete state-space form given by: 

(Xc),+1 = T, + AT * A'.(2/, - 5,e)* (5.3) 

Uk = Xck + A'pJ/2*: 


vhere. 






Figure 5.4: Delta Control Nonlinear Plant 

• AT is the sampling interval. 

r TAS 1 

• j/1 = ^ .is the vector of actual outputs. 

r TASc 1 

• i/c = ti’c j is the vector of commanded outputs. 

• 2/2* the complete vector of state outputs. 

The sampling interval of O.lsec was chosen to be greater than twice the high- 
modal frequency of the uncontrolled plant as shown in Table 4.1. This rate will 
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provide for complete recovery of the analog signal, as per the Nyquist Sampling The- 
orem([Ref. 11]), and is well within the capabilities of the available digital computers. 

Additional modifications included discretizing the differentiators on each state 
and the integrators on each control signal just prior to the controller output. The 
discrete time implementation of these system elements is shown in Figure 5.5 


1^1 

1 ' 

_ 

dscrt_cliff 


1-^1 


dsct1_integ 


Figure 5.5: Discrete Differentiator Sc Integrator 

The final modification results from the fact that the use of a digital control 
induces some time delay in the system. To model this effect a zero-order hold is added 
to the output of the discrete integrator. This block effectively holds the control signal 
for the length of the sampling interval until the next update of the control signal is 
transmitted. Therefore a piecewise continuous control signal is fed to the continuous 
time nonlinear model of the aircraft dynamics. 

The complete discrete controller is detailed in Figure 5.6. It is implemented with 
the nonlinear model in the same fashion as was the differential controller, shown in 
Figure 5.4. 
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Figure 5.6: SiMULiNK Discrete Control Model 

Once the SiMULINK modeled was validated, the discrete controller-nonlinear 
plant system was developed in SystemBuild. The discrete longitudinal and lateral 
controllers are shown in Figures 5.7,5.8. The differentiators and integrators are 
implemented as feedforward and feedback loops respectively. The result is the same 
as for the SiMULiNK implementation shown in Figure 5.5. The blocks incorporating 
the integral gains, also contain the sampling interval. The resulting output executes 
the algebraic expression, 

(sejjt+i = Xk + AT* Kiiy^ - y,)k 
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Figure 5.7: SystemBuild Discrete Longitudinal Control Model 
The blocks entitled 'discJon’ and ’discJat’ complete the discrete controller calcula¬ 
tion by executing the expression, 

Uk = Xck + Kpy2k 

The complete closed-loop model implemented in SystemBuild is shown in Figure 5.9. 

F. MODEL VERIFICATION 

To verify the stability and performance of the closed-loop system, the model 
was subjected to a variety of command and control inputs and external disturbances. 
The state outputs were then evaluated for stability, steady state error, rise time, and 
overshoot. This process was accomplished at each step of the control implementation 
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Figure 5.8; SystemBuild Discrete Lateral Control Model 

process and was the prerequisite for continuation to the next step. Results of the 
final closed-loop time history performance are shown in Appendix C. 

The initialization period of about 10 secs was a result of being unable to set 
initial conditions in the discrete integrators in the SlMULI.XK model. There was 
no such problem when implementing the integrators in SystemBuild. Commanded 
inputs and disturbances were introduced at the following intervals and magnitudes, 

• Commanded TAS increase of 10 ft/sec at time 20 secs. 

• Commanded altitude deacrease of .50 feet at time 40 secs. 

• Commanded heading change of 90 degrees at time 60 secs. 

• Input a headwind of 10 ft/sec at time 80 secs and removed it at time 100 secs. 
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Figure 5.9: SystemBuild Discrete Control Model 

As the plots of Appendix D show, the controller tracks heading (bottom plot of 
Figure D.3), altitude (top plot of Figure D.4), and TAS (bottom plot of Figure 
D.4) with zero steady state error. The feedback system is stable and handles the 
large magnitude heading change by limiting the angle of bank to a maximum of 30 
degrees. This limiting was achieved setting the rate limiter on the heading command 
to +/ —0.2536 rad/sec (approximately 14 degrees/sec). This value was determined by 
examining the dynamics of an aircraft in steady, level, turning flight. The following 
relationships are obtained, 


L * cos4) = mg 
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• L = lift 

• W = mg = weight 

• <j = bank angle 

• R = turn radius 

• V = TAS 



To obtain the maximum turn rate, and thus the bounds for the rate limiter, it is 
necessary only to choose the maximum bank angle and substitute into Equation 5.5. 
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VI. CONCLUSIONS AND 
RECOMMENDATIONS 

A. CONCLUSIONS 

The design process and tools investigated in this thesis allow the designer to 
obtain the optimal controller, model the closed-loop plant/controller system, and 
evaluate the controller’s performance using a computer workstation. The benefits of 
such a process are numerous: 

• Optimal controllers can be designed, evaluated, and modified quickly and easily. 
Simply by adjusting the weighting matrices of the cost function, the designer 
can obtain any desired .system response. 

• Incorporated with a high fidelity model of the aircraft dynamics, the controller 
can be validated throughout a variety of flight regimes. This process proves to 
be cost effective and obtains a level of safety not achievable in actual flight test. 

• Software tools such as MATRIX^ and its graphical design application, Sys- 
temBuild, provide the added benefit of automatically generating the computer 
code required for operational digital flight control systems. The costly and 
time consuming process of generating the thousands of lines of codes required 
by these systems is reduced to development of a block diagram and the click of 
a mouse button. 



B. RECOMMENDATIONS 


Based upon the conclusions and experience gained during this research project, 
the following recommendations are made: 

• The MATR1 X,y / SystemBuild software design package should become an inte¬ 
gral element of the Avionics curriculum at the Naval Postgraduate School. The 
ability to generate real-time C-code makes this combination second to none in 
the area of system simulation and control design. 

• Using the autocode generated from the discrete SystemBuild closed-loop sys¬ 
tem, hardware-in-the- loop testing of the Bluebird controller should be accom¬ 
plished. 

• Incorporate sensor and Kalman Filter models into the SystemBuild block dia¬ 
gram. These models have been developed using SlMULI.NK , [Ref. 4], and must 
be developed for the SystemBuild model to provide a complete model. 

• Chirrently SystemBuild has no block diagram available to model rate limiters 
for the commanded state inputs. These are necessary because the commands 
are usually modeled as step inputs. For large step changes in the commanded 
states, the controller will employ large control deflections to reach the desired 
state. In some instances, this response can be outside the vehicle’s dynamic op¬ 
erating envelope. Rate limiters employed in the SlMl'LINK models have shown 
the capability to limit roll angle, longitudinal acceleration, and vertical accel¬ 
eration to within vehicle operating limits. The same must be done for the 
SystemBuild model. 
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APPENDIX A 


CODED EQUATIONS OF MOTION 


A. STATE_DERIV.M 

function accel = state_deriv(x) 

% calculates the accelerations (angular and linear) due to 
% aero forces; w X v; gravity. 

% Variables brought from workspace: 

% X = combination vector [ contrl inputs, state variables, euler angles] 
% (da, de, dr, dtrt, u, v, w, p, q, r, phi, theta, psi) 

% Variables called from function ’’blue.data” 

% rho = air density 

% b = wing span 

% c = wing cord 

% s = wing area 

% Cfo = Steady state force term 

% Cfu = Stability derivative for control inputs 

% m = airplane mass 

% Ib = inertia tensor matrix (body frame) 

% To = Thrust scale term 
% Pe = Engine postion matrix 

% get the aircraft data 

[ uO,wO,rho,Cfx,Cfo,Cfu,Cfxdot,s,b,c,m,Pe,To,Ib] = blue.data; 


% seperate the combined vector into seperati 
u = [x(l);x(2);x(3)] ; 
dtrt = x(4); 

state = [ x(5); x(6); x(7); x(8); x(9); x(10)] ; 
lambda = [ x(ll); x(12); x(13)] ; 

• ■ [ x(14>, x(15); x(16)] : 


%%%%%% calculate velocity wrt the airmass and form state vector 
%%%%%% that will be used to calculate the aerodynamic forces/moments 

ias = [ state(l); state(2); state(3)] + wind; 

statel = (ias(l)-uO; ias(2); ias(3)-w0; x(8); x(9); x(10)] ; 

%%%%%% calculate total velocity, vt 
vt = (ias(l)A 2 + ias(2)A 2 + ias(3)A 2)A .5; 
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% calculate qbar 

qbar = .5* rho* (vtA 2); 

% calculate Ml 

Ml = diag([ 1/vt, 1/vt, 1/vt, (.5* b)/vt, (.5+ c)/vt, (.5* b)/vt] ); 

% calculate M2 

M2 = diag([ 0, 0, (.5* c)/(vtA 2), 0, 0, 0] ); 

%: calculate Sprime 

Sprime = diag([-1, l,-l,b, c, b] * s); 

% calculate Mu 

Mu = inv([ eye(3)* m,zeros(3);zeros(3),Ib] ); 

9? calculate Tw2b 

alpha = state(3)/vt; 
beta = state(2)/vt: 

T1 - [ cos(alpha), 0, -sin(alpha); 0,1,0; sin(alpha), 0, cos(alpha)] ; 

T2 = [ costbeta), -sin(beta), 0; sin(beta), cosi'beta). 0; 0.0,1] ; 

Tw2b = [ Tl* T2, zeros(3); zeros(3), Tl* T2] ; 

% calculate Chi 

Chi = eye(6) - Mu* T'w2b* qbar* Sprime* Cfxdot* M2; 

% calculate Propulsion matrix 

Prop = [ eye(3); 

0,-Pe(3),Pe(2); 

Pe(3).0,-Pe('l); 

-Pe(2),Pe(l),0] ; 

% calculate gravity vector and rotation matrix 
Rot = [ 1, 0, -sin(lambda(2)); 

0, cos(lambda(l)), cos(lambda(2))* sin(lambda(l)); 

0, -sin(lambda(l)), cos(lambda(2))* cos(lambda(l)j] ; 

Ru2b = [ Rot;zeros(3)‘j ; 
g = [ 0; 0; 32.174] ; 

FgU = m* g: 

% put the components due to gravity; thrust; and control surface deflections 
% together for their contribution to x-dot 

thrust = Prop* To* dtrt; 
gravity = Ru2b* FgU; 

ctrl=qbar* (Tw2b* (Sprime* (Cfo + (Cfu* u)))): 


xdotu=(Mu* (ctrl+thrust+gravity)); 



% calculate rotation matrix 

omegax = [ 0,-state(6),state(5);state(6),0,-state(4);-state(5),state(4),0] ; 
wxlb = (-inv(Ib))* (omegax* Ib); 

Rot = [ -omegax. zeros(3); zeros(3), wxlb] ; 

% rotation component of xdot (w X v) 
xdotrot = Rot* state; 

% state vector feedback component xdot 

xdotcfx =qbar* (Mu* (Tw2b* (Sprime* (Cfx* (Ml* statel))))); 

% add three components of xdot together and premult by inv(Chi) 
xdot= inv(Chi)* (xdotrot-f-xdotcfx-fxdotu); 

% calc accel that a. strap-down IMU would measure 
xdotb=xdot-xdotrot-Ru2b* g; 

% put together for the return vector 
%accd=[ xdotb(l);xdotb(2);xdotb(3);xdot] ; 

return xdot only 

accel=xdot; 


B. EOM2_WIND.C 


Aircraft Equations of Motion; bj' Jim Byerly 
Re\ision 1.1: by Brian T. Foley 


MATRIXx/W'S TM Software V2.4 (C) Copyright 1990 

INTEGRATED SYSTEMS INC., Santa Clara, California 
Unpublished Work. All Rights Reserved Under The U.S. Copyright Act. 

RESTRICTED RIGHTS NOTICE : Use, Reproduction Or Disclosure Is 
Subject To Restrictions Set Forth In The Rights In The 
Technical Data And Computer Clause At 252.227-7013 And The 
Equivalent Provisions In Other Agency's Regulations. 

/* 


This template is provided to let users write their own C 

Code blocks to implement User Code Blocks. See the UCB chapter 

in your SystemBuild manual for further information. 
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iinfo[ 0] ( 0 : Error Occurred : { -10* Message Index -|- Status } 

Return Status : { 0:NormaJ | -LWarning | -2;Error } 


Index Error Message 









0 Provide your own message via MATWR. 

1 Division by 0.0 produces infinity. 

2 Raise 0.0 to a non-positive power. 

3 Both arguments to ATAN2 are zero. 

4 ASIN or ACOS argument out of range. 

5 Natural log of zero or negative number. 

6 Square root of negative number. 

7 Incoming data not in range of table. 

8 Raise negative number to non-integer. 

9 Overflow in y = EXP(u) function. 

H..—..--b 

iinfo 1 =1 : Initialize. (First call in SIM). 

iinfo 2 =1 : Compute the state derivatives. 

iinfo 3 =1 : Compute the outputs. 

iinfo 4 : Reserved. 

iinfo 5 : Reserved. 

iinfo 6 : Reserved. 

iinfo 7 : Reserved. 

iinfo 8 =i : Number of Integer Parameter Values, 
iinfo 9] =r ; Number of Real Parameter Values, 
iinfo 101 =1 : Integration Algorithm (lALGORITHM). 
iinfo[ 11] =1 : Inside Linearization Process 

(Jacobian Computations). 

iinfo] 12] =1 ; Inside SIM Initialization Process (Time=0). 
iinfo] 13] : Update with converged state values. 


Rinfo is a real array of dimension 4 with timing and related 
information for the called routine. User code may not modify 
any values in rinfo. 

: Continuous Discrete Triggered 

rinfo Qi : Current Time Current Time Current Time 

rinfo l] : 0.0 Sample Interval 1.0 

rinfo 2] : 0.0 Initial Time Skew Timing Requirement 

rinfo] 3] : reltol reltol reltol 


The input arguments are defined as follows. Note that initial 
conditions defined in the UCB block form are passed in through 
X, the state vector. The state vector can be modified only on 
the first call to the UCB, when init = 1. LIN UCB‘s are never 
called with init = 1. Hence, the operating point can only be 
changed from usrOl. 

uf nul = input vector dimensioned number of inputs (nu). 

x[ nxj = state vector dimensioned number of states (nx). 

xdotf nx] = first derivative with respect to time of the 
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state vector. If the UCB is in a discrete Super¬ 
block, xdot represents xnext, the value of x at 
the next sample time. 

y[ ny] = output vector, dimensioned number of outputs (ny) 

rpar[ nr] = general vector of real double precision parameter 
(double precision R* 8), to be initialized in 
MATRIXx and passed to the UCB. Nr is requested in 
the SystemBuild Block form of the UCB as ‘dimension 
of RPAR'and passed to the UCB as iinfo[ 9] . 
ipar[ ni] = general vector of integer (1+ 4) parameters to be 

initialized in MATRIXx and passed to the UCB. Ni 
is requested as ‘dimension of IPAR‘in the block 
form of the UCB in SystemBuild, and passed to the 
UCB as iinfo[ 8] . 


I The following two lines of code are not to be modified. 


/ 


int init, state, output, messg; 
double time, tsamp, tskew; 


Enter user variable declarations 


double vbo[ 3] , wbo[ 3] , lambda[ 3] , delu[ 3] , delt, phi, theta, psi; 
double vt, alpha, beta, qbar, mass, s[ 6] , sbar[ 6] [ 6] , m-prinie[ 6] [ 6] ; 
double m.dot.primej 6] [ 6] , mi[ 6] [ 6) , miinv[ 6] [ 6] , chi[ 6] [ 6] ; 
double chiinv[ 6] [ 6] , invnnertia] 3] [ 3] , wcv[ 3] , iwci[ 3] , fvwobf 6] ; 
double wghtf 3] , fgrav[ 6] , deLvw[ 6] , faero[ 6] , fcntl[ 6] , vwnew[ 6] ; 
double slamf 3] [ 3} , euangl 3] , rb2u[ 3] [ 3] , dpos[ 3] ; 
double tmpl[ 6] , tmp3[ 6] [ 6] , tmp5[ 3] , windx, windy, windz; 
double a, b. c; 

double wt = 55.0, rho = .0023769. warea = 22.38. span = 12.42; 
double chord = 1.802, grav = 32.174, TO = 15.0, VO = 73.3, alphO = 0.0 
double cfx[ 6] [ 6] , cfu[ 6] [ 3] , cfx-dot[ 6] [ 6] , cf0[ 6] , inertia] 3] [ 3] ; 


The following six lines of code are not to be modified. 


/ 




= (iinfof 3] 
= rinfo[ 0] 
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tsamp 

= rinfof 1] ; 



tskew 

= rinfo[ 2] ; 



Replace 

the following with us 

er code. 



/* define aircraft data * / 

/* Zeroize Arrays * / 
for (i=0; i( 6; i++) 
for (j=0; j( 6; j++) 

cfx[i][j] =0.0; 

cfx-dot[ i] [ j] = 0.0; 

cf0[ i] = 0.0; 

for (i=0; i( 6; i++) 
for (j=0; j( 3; j++) 
cfu[ i] [ j] = 0.0; 

for (i=0; i( 3; i++) 
for (j=0: j( 3; j++) 
inertia[ i] [ j] = 0.0; 

/+ define elements of arrays * / 

/♦ * * + COLUMNS: U. ALP, BETA. P. Q, R ♦ / 

I* *** ROWS; CD, CY, CL, Cl, Cm, Cn * / 


cfxf 0] 2] = 

0.188; 

cfx 1 1 = 

-0.31; 

cfx ]' 5 = 

0.0973; 

cfx 2 2 = 

4.22; 

cfx 2' 4 = 

3.94: 

cfx 3 1 = 

-0.0597; 

cfx 3 ; 3 ; = 

-0.363; 

cfx 3 5 = 

0.1; 

cfx; 4' ; 2 = 

-1.163; 

cfx 4 ' 4 = 

-11.77; 

cfx' 5 ' 1 = 

0.0487; 

cfx ’ 5 ' ’3 = 

-0.0481; 

cfx[ 5] [ 5 = 

-0.0452; 


/* * * * COLUMNS: da, de, dr ♦ 

I* *** ROWS; CD, CY, CL, Cl, Cm, Cn ♦ / 

cfu[ 0] [ 1] = 0.065; 

cfu 1 2 = 0.0697: 

cfu[ 2] [ 1] = 0.472; 
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cfuf 2 
cfu 3 
cfu 3 
cfu I 4 
cfu 5 
cfu 5 


2] = 0.0147; 

0 = 0.265; 

2 = 0.0028; 

1' = -1.41; 

O' = -0.0347; 

2 = -0.0329; 


U * * * COLUMNS; U, ALP. BETA, P, Q. R ♦ / 
/* + * + ROWS; CD, CY, CL, Cl, Cm, Cn * / 


cfx-dot[ 2] [ 2] = 1.32; 

cfx-dot[ 4J [ 2] = -4.7; 

/* * * * ROWS: CDO, CYO. CLO, CIO, CmO, CnO * / 

cfO[ 0] = 0.03; 
cfO[ 2] = 0.385: 

I* Inertial elements ♦ / 

inertial 0] f 0] = 12.-58; 
inertial Ij I IJ = 13.21; 
inertial 2] [ 2] = 19.99: 


/* Split up input variables * / 


for (i=0; i( 3: i-t-V) 


vboM] 
wbo[ l] _ 
lambda! i 
delu[ ii 

delt = ii[ 12j ; 
windx = u[ 13] ; 
windy = u{ 14] ; 
windz = u[ 1.5| ; 




phi = lambda! 0] ; 
theta = lambda! 1] ; 
psi = lambda! 2] ; 


a = vbol 0] - windx; 
b = vbol Ij - windy 
c = vbof2f- • ’• 
t = sqrt(, 


idz; 

-b b* b -f c 


alpha = asin(vbo! 2] /vt); 
beta = asin(vbo! 1] /vt); 


qbar = .5 * rho ♦ vt* vt; 
mass = wt / grav; 

/* Build S Bar Matrix * / 
s] 0] = -warea; 
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sfl] 
s; 2' 
s 3' 
s 4 
s[ 5 


warea; 

-warea; 

warea * span; 

warea * span; 


for (i=0;i( 6;i++) 


for (j=0; j{ 6; j++) 
sbarf i] f jl = 0.0; 
sbar[ i] [ ij = s[ i] ; 


/* Build m-prime matrix * / 

= l./vt; 

= l./vt; 

= l-/vt; 

= chordV?2''*'n); 

= span / (2.* vt); 
for (i=0;i( 6;i++) 

for (j=0; ]( 6; j++) 
m.primel i] f j] = 0.0; 
m_prime[ i] [ i] = s[ i] ; 

/* Build m-dot .prime matrix * / 

0] = 0.0; 

4 = 0.0; 

2 = chord / (2.* vt* vt); 

3 = 0 . 0 : 

4 ^ = 0 . 0 ; 

5] = 0.0; 


for (i=0;i( 6;i++) 

for (j=0; j( 6; j++) 
m.dot.primel i] [ j] = 0.0; 
m-dot-prime[ i] [ iJ = s[ i] ; 


/* Build inassJnertia matrix * / 
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for (i=0;i( 6;i++) 

j( 

L (i=3; i( 6; i++) 

for 0=3; j( 6; j++} 

mi[ i] 1 j] = inertial i-3] [ j-3] ; 

/* Compute the inverse of the matrix ♦ / 
n=6; 

for (i=0; i( n; i++) 

for 0=0; j{ n; j++) 
ainvfi][j]=mi[i][j]; 

inverse(n); 

for (i=0; i( n; i++) 

for (j=0; i( n: j++) 
mnnv[i] [ j] = amv[ i] [j+n] ; 


/* compute rotation matrices; wind-to-body; inertial-to-body * / 

rot-u2b(phi, theta, psi); 
rot-w2b(alpha, beta); 


!* Compute Chi * / 


rnjnult_666(rfx_dot, m-dot_prime, tmp3); 
m_mult-666(sbar, tmp3, chi); 
m-mult-666(tw2b, chi, tmp3j; 
m-mult_666(rniinv, tmp3, chi); 
for (i=0: i{ 6; i++) 


for (j=0; j( 6; j++) 
chi[ !][]] = -qbar * 


chi[i] [j] ; 


chi[ i] [ i] = 1. + chi[ i] [ i] ; 


/* Compute the inverse of the matrix * / 


n=6; 

for (i=0; i( n; i++) 
inverse(n); 


for (i=0; i( n; i++) 

for (j=0;j(n;j++) 
chiinv[ i] [ j] = ainv 


[i][j+n]; 
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/* Compute w XV*/ 
cross_prod(wbo, vbo, wcv); 

/* Compute inverse of inertia matrix * / 
n=3; 

for (i=0; i{ n; i++) 

for (i=0; j( n; j++) 
ainv[ i] [ j] = inertia[ i] [ j] ; 

inverse(n): 
for (i=0; i( n; i++) 
for (j=0; n; j++) 
invJnertia[ i] [ j] = ainv[ i] [ j+3] ; 

/* Compute I\ (w x Iw) * / 

mjBultJSl(inertia, wbo, iwd); 
cross_prod(wbo, iwci, tmp5); 
m-mult-331(invJnertia. tmp5, iwci); 


for (i—0; i( 3; iH—h) 

fvwobf i] = -wcv[ i] 
fvwob[ i+3] = -iwci[ i] ; 


/* Compute forces due to gravity + / 
wght[ 0] = 0.0: 
wght 1 = 0.0; 
wght[ 2] = wt; 

mjnult_631(tu2b, wght, fgrav); 


/* Compute aero forces and moments * / 


deLvw 

0 

= vbo[ 0] - ■ 

deLvw 

1 

= vbo 1 ; 

deLvw 

2 

= vbo[_2] - 

deLvw 

3 

= wboi 0] ; 

deLvw 

4 

= wbo[ l] ; 

deLvw 

5 

= wbo[ 2] ; 


cos(alphO); 

sin(alphO); 


mjiiultj661(m.prime, deLvw, tmpl); 
m_mult_661(cfx, tmpl, faero); 


mjnultj631(cfu, delu, fcntl); 


for (i=0; i( 6; i++) 

faero[ i] = (faero[ i] + cf0[ i] + fcntl[ i] ) * qbar; 
m_multj661(sbar, faero, tmpl); 
m_multj661(tw2b, tmpl, faero); 
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j* Compute linear (u,v,w) and angular (p,q,r) accelerations + / 

for (i=0; i( 6; i++) 

tmpl[ i] = fgrav[ i] + faero[ i] ; 
t:mpl [ 0] = tmpl[ 0] + TO * delt; 
m_miiltJE)61(miinv, tmpl. vwnew); 

for (i=0; i( 6; i++) 

tmpl[ i] = fvwob[ i] + vwnew[ i] ; 

m_mult_661(chiinv, tmpl, vwnew); 

/+ Compute euler angle derivatives * / 

slamfOl 0= 1.0; 

slam 0 1 = sin(phi) * tan(theta); 

slam 0 2= cos(phi) * tan(theta): 

slam 1 0 = 0.0; 

slam 1 1 = cos(phi); 

slam 1 2 = -sin(phi): 

slam 2 0 = 0.0; 

slam 2 1= sin(phi) / cos(theta); 

slam 2] 2] = cos(phi) / cos(theta); 

m jnull_331(slani, wbo, euang): 

/* Compute position derivatives * / 

for (i=0: i( 3: i-1—t-j 
for (j=0; j( 3: j++} 
rb2u[ i] [ j] = tu2b[ j] f i] ; 

ni_mult-331(rb2u, vbo, dpos); 

/* Set up output vector * / 
for (i=0; i( 6; i++) 

for (1=0; if 3; iT+) 

y[ i+O = euang[ i] ; 
for (i=0: i( 3; i++) 
y[ i+9] = dpos[ i] ; 

} /* End of main program * / 








APPENDIX B 


OPEN-LOOP PLANT 

A. OPEN-LOOP PLANT 


Columns 1 through 7 
-0.0684 0.0000 0.2244 

-0.0000 -0..3877 0.0000 

-0.8621 -0.0000 -4.7.588 
0.0000 -0.1149 -0.0000 

0.0132 -0.0000 -0.2362 

-0.0000 0.0.590 0.0000 

0 0 0 
0 0 0 
0 0 0 
1.0000 0 -0.0000 
-0.0000 1.0000 -0.0000 
0.0000 0 1.0000 
Columns 8 through 12 
-32.1740 0 0 0 0 

0.0000 0 0 0 0 

0.0008 0 0 0 0 

-0.0000 0 0 0 0 

-0.0000 0 0 0 0 

0.0000 0 0 0 0 

0.0000 0 0 0 0 

0 0 0 0 0 

-0.0000 0 0 0 0 

-0.0004 -0,0004 0 0 0 

0.0000 73.3000 0 0 0 

-73.3000 0 0 0 0 


0 0.0021 0 0.0000 
0.0023 0 -72.6109 32.1740 

0 67.9933 0 -0.0002 

4.3390 -0.0000 1.1953 0.0000 

0 -3.8650 0 0.0000 

0.3617 -0.0000 -0.3400 -0.0000 

1.0000 -0.0000 -0.0000 -0.0000 
0 1.0000 -0.0000 -0.0000 
0 0.0000 1.0000 0.0000 
0 0 0 0 
0 0 0 0.0023 

0 0 0 0.0000 



8.7745 


b = 

-5.4349 -0.0000 0 

0.0000 5.8266 0 

-38.7398 -1.2065 0 

-0.0000 0.3949 37.3882 

-26.8913 0.0185 0 

0.0000 -2.9212 -3.0817 

0 0 0 

0 0 0 

0 0 0 

0 0 0 

0 0 0 

0 0 0 


0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 


B. OPEN-LOOP STEP RESPONSE PLOTS 


llie open-loop S 3 'stemBuild model was subjected to a step elevator input of 
approximately 3 degrees a! a time of 31? seconds. TJie following plots detail tfee 
open- loop response of the plant model. The sinusoidal inertial velocitj'. u. .s-bc-vvn in 
the top diagram of Figure B.ldemonslrates the underdamped phugoid mode. .\!so. 
the unstable spiral mode is evident from the linearly ine.reasing plot of heading, h'. 
in the bottom diagram of Figure B.3. 
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APPENDIX C 


CONTROLLER 

CALCULATION /VERIFICATION 

A. BLUE_PLOT.M 

%: This file plots bode diagrams and step responses for the Iqr 
%controller calculated . 

load linl6 

[ alat,blat,alon,blon] =latlonl(a); 

'/o'-X'llliiCtMTipute loagitaiclinal,COTtjolhir'S.Sift 
[ as.bs.cs.ds] =linmod(‘blue-syn. lon‘); 

r=[ 110000 0:0 10000] : 
q=| 1 0:0 1] : 

bl=bs(:.l:2): 

b2=bs(;.3:4): 

cl=as(6:7.;); 

[ k,s,e] =lqr(as.b2,cs‘* q* cs,r); 

damp(e) 

k=2* k; 

kp=k(:,l:5): 

ki=k(:,6:7); 

pause 

%cornputc control bode diagrams 
bode{as-b2* k,b2,k(l,:),[ 0 0] ,1) 
title(‘ele cmd loop') 

bode(as-b2* k,b2,k(2,;),[ 0 0] ,2) 
title(‘thr cmd loop') 

%compute command bode diagrams and step responses. 
step(as-b2* k,bl,cl(l,:),[ 0 0] ,2) 
title(‘step z response') 

bode(as-b2* k,bl,cl(L:),[ 0 0] ,2) 
title('z cmd loop') 



step(as-b2* k.bl xl(2,:),[ 0 0] ,1) 
litle(‘step u response*) 

bocle(as-b2+ k,bl,cl(2,:),[ 0 0] ,1) 
title(‘u cmd loop") 

%%%Gompute lateral controller%%% 

[ as_ lat.bs^ lat,cs_ lat,ds_ lat] =linmod(‘blue- syn_ lat‘); 

r-^ lat=7000; 
q- lat=25; 

bl- lat=bs. lat(:,l); 
b2-lat=bs- lat(;,3); 

cl - lat=as_ lat(l,:); 

[ k- lat,s_ lat,e- lat] =lqr(as_ lat,b2- lat,cs. lat‘* q. lat* cs- lat.r- lat); 

damp(e_ lat) 

k- lat=2* k- lat; 

kp_ lat=k- lat(;,2:6); 

ki. lat=k- lat(:,l); 

compute control bode diagrams 
bode(as_ lat-b2- lat* k_ lat,b2- lat,k- lat,0) 
title(‘aileron cmd loop") 

^compute command bode diagrams and step responses. 
step(as_ lat-b2_ lat* k_ lat.bl. lat,cl_ lat,0) 
title!*step heading response*) 

bode(as_ lat-b2_ lat* k_ lat,bl. lat,cl- lat,0) 

title(‘heading cmd loop*) 

pause 


B. OPTIMAL GAINS 


The .m file shown in the previous section caluclated the optimal gains under the 
assumption of zero coupling between the longitudinal and lateral states. The results 
included a 2x7 matrix for the longitudinal gains and a li7 vector for the lateral gains. 
These gains were then split into proportional and integral error parts for use in the 
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controller gain blocks as detailed in Figure 5.1. The resulting matrices are shown 
below. 

hi.long = 

0.0051 0.0073 
-0.0107 0.0380 
ki-lai = 

0.2689 

kpJong = 

0.0037 0.0095 -0,1632 -1.4319 0.0131 
0.1223 -0.0032 -0,0.3.32 0.0.368 -0.0494 
kpJat = 

0.0077 0.0861 0.2293 0..5062 1.1072 

C. CLOSED-LOOP FREQUENCY AND STEP RESPONSES 
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Amplitude 


step z response 



Time (secs) 

Figure C.3: Step Altitude Command 
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Amplitude 


step u response 



Time (secs) 


Figure C.5: Step Airspeed Command 
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APPENDIX D 


CIOSED-LOOP PERFORMANCE OF 
DISCRETE CONTROLLER ON THE 
NONLINEAR PLANT 


The following plots chronicle the time history of the closed-loop system response 
the commanded inputs and airmass disturbances detailed in Chapter V, Section 
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